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This  paper  reviews  some  recenr.  developments  in  direct  time  integration 
methods  for  nonlinear  structural  dynamics.  The  developments  pertain  to  the 
use  of  linear  multistep  difference  operators  in  conjunction  vich  the  pseudo- 
force approach.  The  paper  is  organized  into  three  main  sections.  An 
introductorv  section  provides  an  overview  of  the  transient  response  analysis 
problem.  A section  on  computational  aspects  deals  with  the  organization  of 
the  numerical  calculations;  this  material  is  largely  based  on  a recent 
detailed  study  of  linear  dynamic  calculations  [1-2).  A section  o;.  integra- 
tion methods  highlights  algorithmic  aspects  that  impact  the  selection  of 
integrator  for  nonlinear  problems  and  discusses  adaptive  analysis  features 
such  as  stepsize  control  and  implicit  matrix  scaling  techniques.  An 
appendix  section  outlines  the  functional  organization  of  modular  "integration 
driving"  software. 
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Section  1 
INTRODUCTION 


1 . 1 The  State  of  Computational  Nonlinear  Mechanics 

Nonlinear  dynamics  is  structural  mechanics' new  frontier.  Although  signifi- 
cant progress  in  the  theory  and  implementation  of  computational  methods  has 
been  made,  nonlinear  dynamics  is  presently  an  unsettled  field  and  is  likely 
to  remain  so  in  the  near  future. 

A decade  ago  over  90%  of  engineering  design  work  and  perhaps  over  half  of 
verification  work  was  carried  out  using  "linear  thinking"  mode.  Of  the 
nonlinear  analysis  fraction,  an  overwhelming  proportion  was  devoted  to  static 
problems  such  as,  for  example,  stability  assessment.  Applications  of  transient 
nonlinear  analysis  were  hampered  by  three  factors : 

(1)  Lack  of  reliable,  computationally-oriented  formulations ; 

(2)  High  computational  costs; 

(3)  Unfamiliarity  of  prospective  users  with  nonlinear  mechanics  in 
general  and  the  use  of  application  programs  in  particular. 

Since  then,  these  roadblocks  have  eroded  in  different  degrees.  Most  dramatic 
progress  has  been  made  in  (1) . Obstacle  (2)  has  weakened  by  gradual  but 
steady  decrease  of  unit  computational  costs  as  well  as  by  improving  efficiency 
of  new  implementations.  The  "education  gap”  (3)  remains,  however,  a formidable 
barrier.  Although  it  is  easy  to  say  that  any  problem  in  structural  mechanics 
can  be  attacked  from  a nonlinear  dynamics  standpoint,  in  practice  a converse 
route  is  followed.  But  it  is  often  difficult  for  the  nonspecialist  to  judge 
in  advance  whether  nonlinearities  and  dynamic  effects  are  of  sufficient 
importance  to  merit  consideration  at  some  stage  of  the  design  or  verification 
cycle.  The  main  purpose  of  this  section  is  to  provide  a quick  overview  of 
structural  problems  that  may  require  a nonlinear  dynamics  treatment. 
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1.2 


Motivation  for  Nonlinear  Analysis 


The  increasing  importance  of  nonlinear  analysis  is  to  a large  extent  due  to 
emphasis  placed  by  manufacturers,  contractors  and  certifying  agencies  on 
realistic  modeling  and  accurate  analysis  of  critical  structural  components. 
This  endeavor  has  prompted  not  only  the  incorporation  of  nonlinear  analysis 
capabilities  into  existing  linear  codes,  but  also  the  development  of  new 
finite-element  and  finite-difference  codes  specifically  designed  for  nonlinear 
analysis . 

The  need  for  nonlinear  analysis  capabilities  was  initially  felt  in  the 
aerospace  industry.  Pressing  demands  for  lightweight  design  led  to  the 
investigation  of  new  structural  configuration  concepts  (e.g.,  thin  shear- 
carrying webs),  new  fabrication  concepts  (e.g.,  honeycomb  panels) , and 
increasing  use  of  brittle  materials  such  as  fiber-reinforced  composites. 
Implementation  of  these  advanced  concepts  prompted  concern  over  realistic 
modeling  of  structural  components  as  regards  numerical  simulation  of  their 
behavior  under  critical  loading  and  environmental  conditions. 

Continuing  progress  in  variational  discretization  and  solution  techniques, 
coupled  with  the  dissemination  of  production-level  software,  has  helped  in 
spreading  interest  in  nonlinear  structural  analysis  to  non-aerospace 
applications  occurring  ir.  automotive,  bio,  civil,  industrial,  nuclear  and 
petroleum  engineering.  Among  those  systems  that  are  often  subject  to  signi- 
ficant nonlinear  effects,  we  can  mention  thin  shells,  pressure  vessels, 
nuclear  reactors,  cable  and  pneumatic  structures,  impacting  structures, 
equipment  mounts,  shock-excited  underwater  and  underarour.d  structures,  and 
metal-forming  machinery. 

1. 3 Static  Versus  Dynamic  Treatment 

As  noted  previously,  an  overwhelming  proportion  of  engineering  applications 
of  nonlinear  structural  analysis  are  performed  in  "static  mode"  despite  the 
undeniable  presence  of  dynamic  effects  in  many  situations,  such  as  snap- 
through  analysis.  This  practice  is  commonly  justified  on  various  grounds. 
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First,  it  can  be  argued  that  many  of  the  problem  characterization  parameters 
such  as  excitation  and  constitutive  behavior  are  presently  poorly  understood 
and  qualitatively  uncertain  to  begin  with.  Second,  the  neglect  of  inertial 
effects  often  leads  to  conservative  results.  Third,  dynamic  analysis  is 
both  more  cumbersome  to  carry  out  and  more  expensive  than  static  analysis. 

As  more  experience  is  gathered,  the  third  argument  weakens.  In  fact,  non- 
linear dynamic  analysis  is  more  easily  implemented  in  a general  purpose 
program  framework  than  nonlinear  static  analysis.  (This  does  not  apply,  of 
course,  to  linear  analysis.)  The  only  practical  solution  procedure  for  the 
former  is  direct  time  integration,  whereas  for  the  latter  there  is  a vast 
array  of  special  techniques  that  must  be  mastered  by  potential  users.  If 
computational  costs  continue  to  decrease  in  relation  to  manpower  costs,  it 
is  not  farfetched  to  presage  that  nonlinear  static  solution  procedures  could 
be  eventually  embedded  into  a dynamic  analysis  driver.  If  this  occurs,  the 
second  argument  becomes  irrelevant.  The  first  argument  (poor  problem 
characterization)  remains  valid,  but  an  answer  to  this  question  lies  primarily 
within  the  scope  of  experimental  mechanics. 

1.4  Sources  of  Nonlinearities 

Problems  in  structural  mechanics  are  formulated  in  terms  of  three  spatial 
fields:  displacements,  strains,  and  stresses.  These  fields  are  connected  by 
the  stress-displacement  and  stress-strain  (constitutive)  relations  as  shown 
in  Figure  1,  which  is  adapted  from  13].  The  boundary  value  problem  is  closed 
by  specification  of  the  boundary  conditions  (B.C.)  on  displacements  and/or 
stresses  (or  stress  resultants).  In  dynamic  analysis,  the  initial  value 
problem  is  completed  by  specification  of  the  initial  displacements  and 
velocities  at  a reference  state. 

Any  of  the  four  two-way  links  shown  in  Figure  1 may  be  nonlinear.  Conse- 
quently, sources  of  nonlinearity  can  be  naturally  classified  into  four 
groups:  force,  material,  geometric,  and  kinematic,  which  correspond  to 

links  1,  2,  3,  and  4,  respectively. 
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Governing  Equations  of  Motion 

The  response  of  the  physical  structure  to  external  excitation  is  approximated 
by  that  of  a discrete  mathematical  model  with  degrees  of  freedom,  which 
are  collect  - in  a state  vector  u.  The  components  of  u are  generalized  dis- 
placement coordinates  that  define  the  spatial  configuration  of  the  discrete 

model  with  respect  to  a reference  state  u for  every  time  t.  (This  reference 

— o 

state  may  be  updated  during  the  solution  of  highly  nonlinear  problems.) 

A contiguous  family  of  conf igurations  u - u(t)  satisfying  initial  and  pres- 
cribed displacement  conditions  is  called  a motion.  The  governing  equations 
of  motion  appropriate  to  the  description  of  a wide  range  of  nonlinear  dynamic 
problems  can  be  presented  as  follows  [ ^ ) : 

Mu  + Du  + <I(u,u)  + JC  u + S(u)  = f_(t)  (1.1)  j 

1 

where  a superposed  dot  denotes  temporal  differentiation,  and 

M,  D,  K mass,  reference-damping,  and  reference-stiffness  matrices, 

Orf 

respectively  (the  "reference’’  concept  is  discussed  in 
Section  1.6)  ; 

C,  £ nonlinear  damping  and  stiffness  operators,  respectively, 

which  generate  state-dependent  corrective  force  vectors; 

Jf(t)  vector  of  applied  (generalized)  forces. 

The  motion  u(t)  corresponding  to  a given  £(t)  will  be  called  the  response  of 
the  system  (1.1).  Essential  B.C.  on  specific  components  of  u(t)  are  assumed 
to  be  embodied  in  (1.1)  for  notational  convenience.  The  effect  of  ncn- 
homogeneous  essential  B.C.  can  be  also  incorporated  in  (1.1)  through  various 
techniques  such  as  the  adjunction  of  penalty  terms  or  the  introduction  of 
Lagrange  multipliers . 


The  generation  of  system  (1.1)  from  the  continuous  field  model  schematized  in 
Figure  1 will  not  be  dealt  with  in  this  paper,  which  is  concerned  with 
techniques  for  solving  it. 


Remark  1 . 


Equation  (1.1)  does  not  cover  the  esoteric  class  of  mass- 


varying  problems  such  as  reentry  vehicle  iynaxnros. 

Remark  2 . Equation  (1.1)  as  written  does  not  model  "rath  dependent" 

problems,  in  which  the  response  is  a functional  or  the  past  motion.  These 
problems  can  be  accounted  for  by  making  the  nonlinear  terms  memory- 
functional  operators. 

Remark  3 . In  explicit  time  integration,  K = 0,  ana  all  stiffness 

information  is  transmitted  in  vector  form. 


1 . 6 The  Pseudo-Force  Approach 

Implicit  direct  time  integration  procedures  for  solving  (1.1)  can  be  based 
on  either  the  tangent  stiffness  or  a pseudo-force  approach.  In  the  former, 
which  is  briefly  discussed  in  Section  2.3,  (1.1)  is  linearized  at  each 
computed  state  and  the  solution  advanced  incrementally.  In  the  latter,  the 
reference  matrices  X and  D are  keDt  fixed  over  mar.v  computational  steDs  anc 
the  nonlinear  terms  S_  and  C are  treated  as  state-dependent  reaction  forces ; 
this  is  the  approach  emphasized  herein. 

The  application  of  the  pseudo-force  approach  to  transient  nonlinear  analysis 
was  first  advocated  by  Stricklin,  Haisler  and  coworkers  i 4 J . Their  strategy 
was  to  use  the  linear  matrices  fC  and  D evaluated  at  the  initial  cor.f iauration 

A/ 

We  can,  however,  do  much  better  from  a global  point  of  view  if  the  selection 
of  X and  D is  not  so  restricted;  in  fact,  such  matrices  need  not  be  identi- 
fiable  with  the  stiffness  and  damping  matrices  evaluated  at  actual  configura- 
tions. This  generalized  pseudo- force  method,  presented  by  Underwood  and 
?ark  [5],  is  designed  to  work  effectively  in  conjunction  with  the  matrix 
scaling  technique  described  in  Section  3. 

The  generalized  stiffness  selection  concept  may  be  illustrated  with  the 
bilinear  spring  shown  in  Figure  2.  In  the  conventional  pseudo-force 
approach,  the  stiffness  would  be  based  on  the  initial  slope  k . In  the 
approach  described  here,  an  average  slope  k can  be  used  to  evenly  distribute 
the  force  residual  magnitude  over  the  expected  operational  range  (note  that 
k is  not  physically  attainable). 


FORCE 


Fig.  2 Reference  Stiffness  Selection  for 

Bilinear  Spring  Element  (after  [5]) 


1-7 


A helpful  guideline  is:  integration  procedures  based  on  the  pseudo-force 
approach  experience  less  difficulties  with  "softening"  systems  (e.q.,  elasto- 
plastic  materials)  than  with  "hardening"  systems  (e.g.,  tensile  structures). 

It  is  therefore  best  to  err  on  the  side  of  overestimating  the  initial  refer- 
ence stiffness.  As  a simple  application  of  this  rule,  consider  the  impact 
analysis  of  two  structures,  and  S^»  which  are  initially  separated.  Then 
the  reference  stiffness  should  be  that  of  S^+S.,  moving  together. 

1 . 7 Classification  of  Nonlinear  Dynamics  Problems 

A simple  classification  scheme  can  be  based  on  the  effect  of  the  spectral 
characteristics  of  the  excitation  on  the  overall  structural  response: 

Structural  Dynamics  (SD)  Problems:  The  response  is  controlled  by 

a relatively  small  number  of  low-frequency  modes. 

Wave  Propagation  (WP)  Problems:  The  contribution  of  intermediate 
and  high-frequency  structural  modes  to  the  response  remains 
important  throughout  the  time  span  of  interest. 

In  the  above,  the  terms  "frequencies"  and  "modes"  refer  to  the  eigenspectrum 
of  the  linearization  of  the  left-hand  side  of  ( -1) ; because  of  nonlinearities 
this  spectrum  can  be  expected  to  vary  (usually  mildly)  with  the  state.  The 
qualifiers  low,  intermediate,  and  high  refer  to  frequencies  whose  associated 
wavelengths  are  much  larger  than,  of  the  order  of,  and  much  smaller  than  the 
characteristic  acoustic  wavelengths,  respectively. 

In  practice,  a combination  of  the  preceding  types  is  often  encountered: 

Shock-excited  Problems:  An  initial  high-frequency  transient 
gradually  decays  to  a steady  state  or  free  vibration  regime  (step 
waves,  earthquakes,  single  impact,  separation  phenomena,  accidents). 

Multishock-excited  Programs:  Multiple  high-frequency  transients 
occur  with  some  regularity,  possibly  leading  to  gradual  collapse 
or  fatigue  phenomena  (reflecting  shock  waves,  repeated  impacts, 
periodic  deployment) . 
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Dynamic  Instability:  Relatively  smooth,  bounded  response 
suddenly  "takes  off"  as  critical  conditions  are  attained 
(follower  forces,  flutter,  parametric  excitation,  phase 
changes) . 

Generally  speaking,  implicit  integration  methods  are  most  effective  for  SD 
problems  (or  phases)  while  explicit  integration  methods  are  at  their  best 
for  tracing  wp  problems  (or  phases)  . 

1.8  Organization  of  the  Paper 

Direct  time  integration  techniques  for  (1.1)  can  be  studied  from  two 
standpoints : 

1.  Algorithmic : Study  of  mathematical  properties  such  as  stability 

and  accuracy  of  the  integration  formulas. 

2.  Computational : Study  of  properties  directly  correlated  with  the 

implementation  of  the  integration  procedure  on  a computer. 

An  array  of  subjects  pertaining  to  the  latter  viewpoint  is  presented  in 
Figure  3 in  logical  flowchart  format.  The  basic  formulation  steps  leading 
to  the  "advancing  cycle"  box  (central  column  in  Figure  3)  are  covered  in 
Section  2. 

Section  3 reviews  algorithmic  aspects  that  impact  the  selection  of  integra- 
tion formulas  for  nonlinear  dynamics,  and  discusses  features  labeled  as 
"advanced"  in  Figure  3,  namely  accuracy  control  via  automatic  stepsiza 
selection,  and  the  matrix  scaling  technique. 
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Section  2 

COMPUTATIONAL  ASPECTS 
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2.1  Reduction  to  a First-Order  System 


Equation  (1.1)  can  be  reduced  to  a first  order  system  by  introducing  the 
auxiliary  vector  v ( 6 ) : 

v=AMu+Bu  (2.1) 

. — — 

v * AMu  + Bu 
— — 

( 2.2 ) 


where  A and  jJ  are  (for  the  moment)  arbitrary  n^  by  n^  matrices  with  A non- 
singular. These  matrices  are  independent  of  time  and  state,  although  they 
can  be  functions  of  the  integration  stepsize.  The  resulting  first-order 
system  (2.1  - 2.2)  can  be  written 


(2.3) 


where  ^ and  jO  denote  the  identity  and  null  matrix,  respectively,  and  Q 
is  the  effective  force 

2 = (2.4) 

Note  that  the  nonlinear  terms  C and  S have  been  transferred  to  the  right- 
hand  side  as  "pseudo-forces."  Equation  (2.3)  fits  within  the  standard  O.D.E. 
format 

v = F(v,t)  (2.5) 
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2 . 2 Time  Discretization  by  LMS  Formulas 


(2.6) 


We  consider  the  numerical  integration  of  the  first-order  system  (2.3)  using 
m-step , one-derivative,  linear  multistep  (LMS)  difference  operators  to  dis- 
cretize the  time  variation  of  u and  v,  A detailed  treatment  of  these  methods 
can  be  found  elsewhere,  e.g.,[7-9).  For  a constant  stepsize  h,  LMS  operators 
applied  to  Eq.(2.J)  at  the  n-th  time  station  t can  be  written 
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where  o^,  B^,  and  5^  are  scalars  associated  with  specific  finite 

difference  operators,  and  u,  v\  , ...  stand  for  the  vectors  u(t  ) , v(t  ) , ...  , 

K K K jv 

computed  through  the  discretization  procedure  (2.7).  The  above  difference 
expressions  may  be  scaled  by  arbitrary  factors.  In  the  sequel  they  will  be 
normalized  by  requiring  that 


a = Y = 1 (2.8) 

o o 

If  Sq  is  nonzero  (zero),  the  operator  (2.7a)  is  said  to  be  implicit  (explicit) , 
and  similarly  for  (2.7b). 


With  a view  to  subsequent  manipulations,  it  is  convenient  to  introduce  a 
more  compact  notation  which,  nevertheless,  separates  explicitly  the  current 
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fig*,  .^tegggga 


state  variables  u , 

—n 


solutions : 


u , 
— n 


v , v from  historical  terms  involving  past 
— n — n 


u 

— n 

+ 

u 

a 

— n 

= 

hs 

u 

— n 

+ 

h 

bu 

— n 

V 

— n 

+ 

V 

c 

— n 

= 

h5 

V 

— n 

+ 

h 

— n 

where  the  "generalized  stepsizes"  h , h^  are  defined  by 
h„  = h 6 h.  = h 5 

PO  6 O 

and  vectors  <a,  b,  ...  denote  linear  combinations  of  past  solutions: 

al  6 i Y i 

a 8 Y 6 
m m m m 


Z 

, z 

2 

, 2 

a 

b 

c 

d 

2 _ 

2 

— n 

—n 

— n 

— n 

— n-1 

— n-m 

(2.9) 


(2.10) 


(2.11) 


In  Eq.  (2.11),  z is  an  arbitrary  vector  symbol,  which  may  stand  for  a,  u,  .. 

etc.  The  current  state  solution  (u  ,v  ) can  be  expressed  from  (2.9)  in  the 

— n -n 

fc  l 


/ 

t \ 

1 

• 

, u 

U 

u 

h 

—n 

J 6 

— n 

— n 

= \ 

+ 

I 

• 

, V 

V 

l 1 

V 

h 

— n 

\ « 

— n 

—n 

(2.12) 


where 


hU 

-ti 

= 

h 

bU 

- 

u 

a 

— n 

hV 
— n 

- 

h 

.V 

a 

— n 

- 

V 

c 

—n 

will  be  called,  following  Jensen  [6],  the  historical  vectors  pertaining  to 
u and  v,  respectively. 
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2.3 


Implicit  Integration 


We  consider  first  implicit  LdS  operators,  for  which  both  and  are 

nonzero,  whence  h„  j*  0 and  h,  ? 0.  Elimination  of  u and  v from  (2.3)  and 
8 5 — n — r. 

(2.13)  then  produces  the  nonlinear  algebraic  system  of  order  2 n^: 


A M + h.  B 
' — * 8 — 


-h.  i 
8 ~ 


AD-B+h,  AK  nl 

~ ~ ~ g ~ <■- 


v 
— n 


u 

A M h 
— n 

(A  D - B)  hu  + nhV  + hQ  A Q 
— ~ ~ — n -n  8 ~ 


(2.14) 


where  n = h_/hr  = g /6  , and  Q = f - C - S . Finally,  elimination  of  v 
8 5 o o -n  — n — n 

from  (2.14)  yields  the  r.f-th  order  nonlinear  algebraic  system 


•n 


E u 


= 


(2.15) 


where 


E = M + h,  D + h.  h,  K + (h.  - hJ  G 

~ ~ 5 ~ 8 5 ~ B 5 ~ 

A.  ■ [»  + V2.-S’]  £ * hB  * ' 1 £ * h8  h5  2* 


(2.16) 


a"1  b 


The  selection  of  a common  LMS  operator  for  both  u and  v is  natural  if  the 
analyst  wishes  to  keep  similar  discretization  accuracy  on  both  vectors.  In 
this  case,  Eqs.  (2.9)  become 


u 

— n 

+ 

u 

a 

= 

h3 

u 

— n 

+ 

h 

bu 
— n 

V 

— n 

+ 

V 

a 

— n 

= 

h8 

V 

—n 

+ 

h 

bV 
— n 

(2.17) 


and  the  components  of  (2.15)  simplify  slightly  to 


l - fl  * »e£  * »8  £ 


- [a  +h8  <£-£>]  * hs  i: 1 £ * be2  A, 
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The  coefficient  matrix  .E  becomes  independent  of  the  choice  of  A and  B in 
this  case,  as  noted  by  Jensen  (6  (. 

Remark  1 . The  generality  of  the  discretization  (2.7),  in  which  different 

LMS  operators  may  be  used  for  u and  v,  has  two  important  advantages.  First, 
the  use  of  different  approximations  for  u and  v is  allowed,  should  that  prove 
desirable  on  account  of  the  physical  meaning  of  v resulting  from  specific 
choices  of  A and  £ in  (2.2).  Second,  special  two-derivative  integrators  devised 
for  treating  the  second  order  system  (1.1)  directly,  such  as  the  widely  used 
Newmark  [10),  Houbolt  [11) , and  Wilson  [12 | operators,  can  be  presented  as 
special  cases  of  a slight  extension  of  (2.7),  in  which  a historical  u-term 
is  appended  to  (2.7b),  as  shown  i:i  Reference  (lj.  It  follows  that  a separate 
study  of  the  computer  implementation  of  these  methods  is  not  required. 

Remark  2 . If  h^  ? h^ , the  choice  of  A and  should  be  such  that  G = ^ ^ B_ 

is  symmetric  and  maintains  the  sparseness  characteristics  of  the  structural 
matrices.  (All  of  the  choices  discussed  in  Section  2.5  comply  with  these 
restrictions. ) 

Remark  3.  If  two  different  A-stable  operators  are  used  for  u and  v,  and 

B is  nonzero,  the  resulting  system  (2.14)  is  not  necessarily  A-stable.  Sta- 
bility conditions  have  to  be  checked  out  in  advance  to  avoid  surprises. 

2.4  Computational  Paths  for  Implicit  Integration 

Implementations  of  the  implicit  integration  procedure  outlined  in  the  pre- 
ceding section  may  differ  in  two  respects:  (a)  the  selection  of  weighting 
matrices  A and  B_  in  (2.1)  , and  (b)  the  manner  in  which  the  resulting 
auxiliary  vector  v and  its  temporal  derivative  v are  computed  at  each  step. 

We  deal  with  the  latter  topic  first. 

Three  basic  computational  paths,  labeled  (0),  (1),  and  (2),  may  be  followed 
in  advancing  the  discrete  solution  over  a typical  time  step.  The  associated 
sequences  of  calculations  are  flow  charted  in  Figure  4. 

The  path  identification  index,  (0-2) , gives  an  idea  of  the  number  of  backward 
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— — STARTING  PROCEDURE: 

INITIAL  CONDITIONS,  COMPUTE  u. , u. 
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difference  operations  performed  in  the  determination  of  v and  v in  each 

- n — n 

solution  advancing  cycle.  It  is  shown  in  Reference  |2]  that  this  index 
plays  an  important  role  in  the  error  propagation  behavior  of  particular 
implementations . 

Path  (0)  is  consistent  with  the  first-order  difference  system  (2.14),  i.e., 

the  computed  vectors  satisfy  this  system  exactly  if  computational  errors 

are  neglected  and  the  solution  of  the  nonlinear  system  (2.15)  is  iterated  to 

convergence.  (The  solution  error  can  actually  be  regarded  as  a component  of 

the  ’’local"  computational  error  [ 2 J . } The  original  O.D.E.  system  (1.1)  is 

also  satisfied.  On  the  other  hand,  the  computed  u , u and  v do  not,  in 

— n — n — n 

general,  verify  the  auxiliary  vector  definition  (2.1)  , which  is  satisfied 
only  in  the  limit  h -*■  0. 

A slight  variant  of  path  (0),  labeled  (O’),  is  also  shown  in  Figure  4.  In 
this  variant  the  velocity  vector  u , is  recomputed  so  that  (2.1)  is  verified 

“Tl-1 

at  past  stations.  This  "delayed  correction"  is  not  only  expensive  for  general 
A and  B but  worthless  as  regards  improvement  of  error  propagation  character- 
istics [1,2] . It  is  included  in  Figure  4,  however,  because  it  occurs  naturally 
in  the  conventional  choice  v = u discussed  in  Section  2.5. 

Path  (1)  satisfies  both  differential  expressions  (1.1)  and  (2.1)  at  the 
expense  of  (2.14b).  Finally,  path  (2)  enforces  the  v definition  (2.1)  and 
the  difference  system  (2.14)  at  the  expense  of  (1.1). 

2 . 5 Selection  of  Auxiliary  Vector 

The  computational  effort  per  time  step  (and,  to  a lesser  extent,  the  storage 
requirements)  can  be  significantly  reduced  for  certain  choices  of  A and  B 
because  some  of  the  computational  steps  displayed  in  Figure  4 can  be  either 
simplified  or  entirely  bypassed.  The  following  two  selections  were  studied 
in  detail  in  Reference  [1]  for  linear  dynamic  problems: 

u (A  = M-1  , B = 0)  (2.19) 

Mu  + £ H <£  = £ ' £ = £)  (2.20) 
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Equations  (2.19)  and  (2.20)  were  labeled  the  conventional  (C)  and  Jer.senvs 
(J)  formulations,  respectively.  When  these  two  choices  are  combined  with 
the  four  computational  paths  of  Figure  4,  a total  of  eight  formulations 
of  the  advancing  step  result.  The  comparative  study  carried  out  in 
Reference  (:]  identified  the  following  combinations  as  preferable: 

(Cl)  Minimizes  computational  effort  for  undamped  systems  when  used 
in  conjunction  with  backward-difference  operators  and  requires 
less  vector  storage  than  (JO) , but  does  not  possess  favorable 
error  propagation  control  properties. 

(JO)  Minimizes  computational  effort  per  step  for  generally  damped  systems 
and  has  excellent  error  propagation  control  properties,  but  requires 
additional  vector  storage. 

In  both  these  formulations,  the  presence  of  an  ill-conditioned  or  singular 
mass  matrix  M does  not  cause  any  particular  difficulty.  There  are  also  no 
problems  associated  with  the  starting  procedure.  Similar  conclusions  can  be 
extended  to  transient  nonlinear  analysis  handled  with  the  pseudo-force 
approach.  The  optional  corrective  iteration  cycle  shown  in  Figure  4 does  not 
affect  computational  effort  rankings  provided  invariant  calculations  are 
moved  out  of  the  loop.  The  computational  sequences  associated  with  the  (Cl) 
and  (JO)  formulations  are  shown  in  Tables  1 and  2,  respectively. 

Additional  formulations  of  interest  in  the  nonlinear  case  are  obtained  if 
matrices  A and  _B  are  allowed  to  be  stepsize-dependent . For  example,  the 
selection 


v = Mu  + (D  + h.  K)  u (2.21) 

— ~ — ~ S ~ — 

in  conjunction  with  path  (0)  leads  to  a variant  of  the  (JO)  formulation 
that  will  be  denoted  by  (JOK) . The  associated  computational  sequence  is 
shown  in  Table  3.  Its  more  interesting  feature  is  the  relocation  of  steps 
(a)  and  (b) , which  is  believed  to  offer  some  advantages  in  array  manage- 
ment if  an  iteration  loop  is  applied. 


2-8 


Step 


Calculation 


(2)  If  iteration  is  used,  calculation  of  the  invariant  portion  of 
g should  be  moved  outside  the  loop. 


Table  2 


r 

> 

t 


COMPUTATIONAL  SEQUENCE  ASSOCIATED  WITH  THE  (JO)  FORMULATION 


Step 

Calculation 

a 

V 

3 

Q , - K u , 

-n-l 

%-l n-l 

.v  , 

b 

V 

=S 

n , + h,  v 

-n-l 

— n-l  c n—  1 

V 

, v 

c 

h 

=s 

h d - c 

— n 

— n — n 

d 

hu 

. . u u 

h b - a 

— n 

— n — n 

(1) 

e 

E 

M + h„  D + ilr  h K. 

~ S ~ S 

f 

Predict 

or  re-evaluate  C,S 

-n  -n 

1 

1 

9 

2n 

= 

M hU  + hQ  hV  + hr  h„  Q 
n S -n  5 3-% 

1 

1 

-1 

J 

h 

u 
— n 

£ 2n 

1 

(2) 

, _u,  „ Optional  Iteration 

i 

u 

— 

u - h )/ha  

— 

— n 

-n  -n  S 

MOTES 


(1)  , (2) , see  Table  1 . 


Table  3 


COMPUTATIONAL  SEQU £NCE  ASSOCIATED  WITH  THE  (JOK)  FORMULATION 


i 


Step 


, V 

V 

v 

c 

h 

= h 

d - 

c 

— n 

— n 

— n 

, u 

, u. 

a 

d 

h 

= h 

b - 

a 

— n 

ij 

— n 

~n 

(1) 

2 

e 

E 

- M 

D *.•  ho  K 

Calculation 


Predict  or  re-evaluate  C , S 


a 

V 

— n 

— 

2n  " 

K h 
n 

g 

s 

M hu 
n 

h 

u 

— n 

= 

e'l 

2n 

i 

u 

— n 

= 

(u  - hu)/h, 
— n — n 3 

b 

V 

— n 

=5 

hv  + 
~n 

h.  v 

5 — n 

Outional  Iteration 


(2) 


NOTES:  (1),  (2)  see  Table  1. 
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Similarly,  the  use  of  the  auxiliary  vector 


v = M (u  - h3  u ) = M hU  (2.22) 

— ~ — n S — n ~ — n 

in  conjunction  with  path  (2)  allows  a similar  relocation  of  steps  (a)  through 
(b) ; its  practical  value  has  not  been  assessed,  however. 

Table  4 provides  a "snapshot"  view  of  vector  structures  involved  in  the  pro- 
gramming of  the  (JOK)  sequence.  The  flow  of  computations  required  to  advance 
the  solution  over  the  n-th  step  is  indicated  with  arrowed  paths.  (Similar 
diagrams  for  the  (Cl)  and  (JO)  sequences  can  be  found  in  |1|.) 

2 . 6 Primary  and  Secondary  Computational  Variables 

In  the  formulations  discussed  in  previous  sections,  the  displacement  vector 

un>  which  is  obtained  from  solving  (2.15),  is  called  the  primary  computational 

variable;  the  other  three  state  vectors  (u  , v , v ) are  ;econdary  or  derived. 

— n - n — n — ■ 1 — 

We  could  have  solved  (2.3)  and  (2.12)  for  another  state  vector.  For 

example,  the  counterpart  of  (2.15)  for  the  velocities  is 

= -(h6K  + G)h^  + A'1  h^  + h52n  (2.23) 

The  displacements  u would  then  be  computed  from  (2.12a);  this  sequence 

_n 

effectively  reverses  steps  (h)  and  (i)  in  Figure  4 . (Corresponding  systems 

for  v and  v are  considerably  messier  for  general  A and  B.) 

— n — n ~ ~ 

If  exact  arithmetic  were  used  throughout,  the  selection  of  primary  variable 
would  be  irrelevant.  The  choice  does  matter,  however,  if  computational  error 
propagation  is  considered.  Generally  speaking,  the  state  variable  of  main 
interest  in  the  analysis  should  be  selected  as  primary.  For  structural 
dynamics  problems,  displacements  are  a natural  choice.  In  fluid-structure 
interaction  problems,  however,  structural  velocities  are  of  utmost  concern  and 
have  been  used  as  primary  computational  variables  in  the  implementation  of 
staggered  solution  procedures  ( 13 ] . 


Table  4 

VECTOR  STORAGE  ARRANGEMENT  FOR  (JQK)  FORMULATION 


Multistep  Historical 

Previous  step 

Current  Step 

Information  (missing  if  m=l) 

Information 

Information 

NOTES : 


If  a pair  of  backward  difference  operators  (b  = d = £)  is  used, 
information  enclosed  in  brackets  is  not  required  for  advancing 
the  computations,  and  those  computational  steps  identified  by 
broken-line  paths  may  be  omitted. 

If  iterations  are  applied,  the  arrangement  is  slightly  more 
complex,  since  invariant  vector  portions  must  be  removed  from 
the  iteration  loop. 
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Explicit  Integration 


If  two  explicit  integration  procedures  (6  =5  =0)  are  selected  in  (2.1) , 

o o 

the  computational  sequence  becomes 


u 

— n 

= hu 
— n 

V 

— n 

= hV 
— n 

u 

n 

- M-1  A-1 

(v  - B u ) 

— n ~ — n 

V 

— n 

■ A (^  - 

Ku  - Du)+Bu 

~ — n ~ — n ~ — n 

Economic  execution  of  this  sequence  demands  the  presence  of  a diagonal  (and 
nonsingular)  mass  matrix  M.  In  addition,  the  conventional  formulation 
v - u,  v = ii  JL  = £)  simplifies  the  last  two  steps.  The  third  one 

becomes  trivial,  and  the  only  significant  matrix-vector  operations  occur  in 
the  last  step  (acceleration  calculations) . 

The  "full-step"  sequence  (2.24)  should  be  avoided  in  practice,  however, 
because  of  its  poor  error  propagation  characteristics . The  so-called  "half- 
step" or  "half-station"  formulation  is  much  preferable: 


u 

— n 

= h 
— n 

v , 
—n+h 

= u . 
—n+h 

* h*1 , 

(2.25) 

V 

— n 

- u = 

—n 

M-1  (Q  - K u -Du) 

~ ~ — n ^ — n 

The  summed  form  [7]  of  the  central  difference  method  can  be  presented  in  this 
form  [14].  Difficulties  arise,  however,  if  the  system  is  damped  because  the 
velocity  term  in  (2.25c)  is  not  readily  available  in  this  implementation. 
Reference  [ 1 5 J discusses  various  ways  of  circumventing  this  problem. 
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The  Tangent  Stiffness  Approach 


We  close  Section  2 with  a brief  review  of  the  tangent  stiffness  approach  to 

nonlinear  structural  mechanics.  Assume  that  the  computations  have  proceeded 

up  to  t - t . Introduce  the  increments  to  the  next  time  station: 
n 

def 

Az  = z - z (2.26) 

— n — n+1  — n 

where  z_  stands  for  any  variable  vector  such  as  u,  u,  ...,  etc.  The  effective 
force  increment  is  linearized  as  follows: 


*2n 


Af  - AC  - AS 
— n — n — n 


= Af  - (£  + S ) Au  - C.  Au  + 0( | |Au| I 2)  (2.27) 

— n '•‘u  ~u  — — 11  11 

in  which  C , C.  and  S denote  the  matrix  Jacobians  of  C and  S with  respect  to 
~u  '-a  ~a  — 

u and  u,  evaluated  at  t^.  The  tangent  stiffness  equations  can  be  derived  as 
follows.  Insert  a A in  front  of  every  variable  state  vector  in  (2.3)  and 
(2.12),  replace  A£^  in  (2.3)  by  the  linear  portion  of  (2.27),  eliminate  Au^, 

Av  and  Av  as  in  Section  2.4,  and  finally  divide  through  by  h„h..  The  result 
is  the  "quasi-static"  incremental  counterpart  of  (2.15)  (see  also  [3]): 


K*  Au  = H + Af 
~ — n — n — n 


(2.28) 


where 


it 

K 

= K + C 

+ S 

* h’s1 

(D  + c-)  + (h71-hr1) 

~ o B 

-1 

1 £ + (h3h6)  M 

H 

— n 

- ['Vo1 

I’1  M 

(D-G)] 

*u  —1  — 1 

h + h.  A 

1 -n  o ~ 

hV 
— n 

-n 

= hU  - 

“n+1 

u + 
~n 

h.  u 

S -n 

hv 

“n 

= hV  - 
“n+1 

v + 
-n 

h,  £ 
o -n 

(2.29) 


In  (2.28),  K is  the  tangent  (dynamical)  stiffness  and  is  a historical  force 
vector. 
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la  the  "purely  incremental"  application  of  this  approach,  Eq.  (2,28)  is 

applied  at  each  step  to  advance  the  solution  bv  Au  and  no  correction  is 

- n 

applied  at  t ^ to  account  for  the  truncation  error  in  (2.27).  in  large  two- 

and  three-dimensional  problems,  this  naive  procedure  incurs  in  enormous  com- 

* 

putational  expense  in  forming  and  factoring  K at  each  step  while  achieving 

only  first-order  accuracy  (regardless  of  the  accuracy  of  the  integration 

operator) . A better  strategy  is  to  apply  one  or  more  corrections  at  t , . 

n+l 

This  helps  in  balancing  the  accuracy  of  the  integration  formula  with  that  of 
the  nonlinear  force  calculations. 

* 

If  K is  maintained  constant  over  many  steps  and  recomputed  only  as  needed 
because  of  strong  nonlinearities,  a restricted,  "tangent-updating"  version  of 
the  pseudo-force  method  results.  As  noted  in  Section  1.6,  a more  general 
strategy  is  emphasized  in  this  paper.  By  judicious  selection  of  the  reference 
"secant"  matrices  in  (1.1) , highly  nonlinear  problems  have  been  processed 
without  a single  recalculation  of  the  left-hand  side  matrix  E in  (2.15) . 
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Section  3 


INTEGRATION  METHODS 


3.1  General 


Three  major  aspects  of  integration  methods  for  nonlinear  structural  problems 
are  covered  in  this  section*  stability  and  accuracy  characteristics  of 
integration  formulas,  stepsize  changing  techniques,  and  automatic  stepsize 
selection  strategies.  Throughout  this  section  we  use  the  term  "method"  to 
represent  the  synthesis  of  these  three  aspects. 

The  question  of  stability  of  implicit  formulas  for  nonlinear  problems  is  far 
from  settled  and  will  undoubtedly  continue  to  receive  much  attention.  A 
review  of  some  promising  results  is  given  in  Section  3.2 

The  bulk  of  the  literature  on  integration  methods  is  concerned  only  with 
algorithmic  properties  (stability  and  accuracy)  of  fixed-step  integration 
formulas.  In  practice,  these  properties  are  seldom  preserved  when  the  step- 
size  varies.  Moreover,  in  implicit  methods,  the  dynamic  coeffioient  matrix 
has  to  be  formed  and  factored  whenever  the  stepsize  is  changed.  We  use  here 
the  term  "step  changing  technique"  to  denote  the  set  of  rules  and  procedures 
for  minimizing  the  effect  of  stepsize  changes  on  stability  and  accuracy  of  the 
fixed-step  formula,  and  for  reducing  the  computational  effort  involved  in  each 
stepsize  change.  This  is  dealt  with  in  Section  3.3. 

Automatic  stepsize  selection  strategies  are  the  most  talked  about  but  least 
advanced  of  the  three  aspects.  Two  crucial  questions  in  computational  struc- 
tural dynamics  are:  How  can  ar.  explicit  solution  process  be  kept  stable  when  low 
accuracy  is  requested,  and  how  can  an  apparently  stable  solution  computed  by 
implicit  methods  be  guaranteed  accuracy?  In  Section  3.4,  we  present  a pro- 
mising step  selection  strategy  for  explicit  methods  which  has  been  shown  to 
be  both  robust  and  reliable,  and  which  is  aimed  to  answer  the  first  question. 
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For  implicit  methods  we  also  list  promising  error  control  criteria  and  give 
an  appraisal  based  on  somewhat  limited  experience. 

In  Section  3.5  we  mention  some  of  the  new  approaches  aimed  to  solving 
effectively  special  classes  of  problems.  Finally,  we  list  some  of  the  pro- 
mising areas  for  future  research  in  Section  3.6. 

3.2  Stability  and  Accuracy  of  Multistep  Methods 

Let  us  consider  the  homogeneous,  undamped  nonlinear  case  with  identity  mass 
matrix 


u + £ (u)  =0  (3.1) 

Following  Prothero  and  Robinson  | 16),  we  approximate  (3.1)  in  the  neighborhood 
of  the  numerical  solution  u * £ as 

..  n 2 

u sss  ~2<y)  ~ (w  ) (u  - 2)  (3.2) 

N 2 3 

(<u  ) - 2(u)  (3.3) 

The  (Cl)  implementation  of  (3.2)  using  (2.16),  i.e.,  the  same  integration 
formula  for  u and  u,  yields  after  some  manipulation 
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The  stability  of  multistep  formulas  for  model  nonlinear  equations  (3.1) 
can  be  assessed  from  the  characteristic  equation 


TV-*  a. 

/ > ~i 


0 ; X . « 1 

i 


and  the  local  error  is  obtained  from  (3.4)  by  taking  {e  ,*  0;  j=l,..,  k)  . 

n-3 
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We  now  deduce  some  known  results  from  (3.4  - 3.6)  when  Equation  (3.1)  reduces 

to  a scalar  problem.  For  the  trapezoidal  rule  (a  = -a.  = 1,  3 = 5,  = 0.5), 

o 1 o 1 


Equation  (3.5)  becomes 


(A  - 1)Z  + (X+i)  (a,2  A + u2  ) =0 

4 \ n n-1 1 


from  which  one  obtains  the  local  stability  condition  of  the  trapezoidal  rule 


w > w . 
n — n-1 


The  stability  condition  (3.8)  implies  that  the  trapezoidal  rule  is  locally  un- 
stable for  spring-softening  in  the  sense  that  the  solution  perturbation  grows 
(see  Fig.  5) , viz. , 
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Fig.  5 Local  Error  Amplification  Factor  of  Trapezoidal 
Rule  for  a Moilul  Nonlinear  Problem 
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This  undesirable  property  was  pointed  out  by  Gourlay  |17|  for  first-order 
parabolic  equations  and  for  second-order  hyperbolic  equations  by  Park  (18,19]. 
in  which  a pathological  case  was  presented.  For  bilinear-spring  models, 
several  pathological  examples  were  presented  by  Hughes  | 20  | ,•  detailed  condi- 
tions for  unstable  solutions  were  recently  derived  by  Westermo  | 21 ] by  the 
phase-plane  technique. 


The  local  error  is  easily  computed  from  (3.16) 


from  which  two  asymptotic  cases  are  obtained 
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(3.11) 


The  asymptotic  expressions  (3.11)  show  that  the  local  error  can  be  made  small 
for  the  two  extreme  frequency  components.  However,  nothing  can  be  said  from 
(3.10)  as  regards  the  magnitude  of  the  local  error  corresponding  to  intermediate 
frequency  components.  In  a practical  context,  this  implies  that  even  though  j 

the  stability  condition  (3.8)  is  violated  for  the  two  extreme  frequency  com-  \ 

ponents,  the  absolute  error  magnitude  will  most  likely  remain  small  and  may 
not  be  detected.  This  suggests  that  the  instability  of  the  numerical  solution 
b>  the  trapezoidal  rule  becomes  noticeable  when  it  is  caused  by  intermediate 
frequency  components  as  explained  in  i 22 i . 
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The  source  of  the  instability  of  the  trapezoidal  rule  for  certain  nonlinear 
problems  was  attributed  to  the  use  of  historical  derivatives  in  [IS] . 

The  effect  of  historical  derivfatives  can  be  eliminated  by  adopting  one-leg 
formulas,  which  are  obtained  by  selecting 
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(3.11) 


Notice  that  backward  difference  formulas  are  a special  class  of  the  one-leg 
formulas.  Let  us  consider  the  midpoint  implicit  formula,  which  is  the  one- 
leg  adaptation  of  the  trapezoidal  rule,  viz., 
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The  characteristic  equation  for  the  midpoint  rule  (3.12)  is 
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(3.13) 


which  indicates  stability  regardless  of  the  possible  variation  of  The 

local  error  is  obtained  in  the  form 
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(3.15) 


The  comparison  of  the  two  asymptotic  error  expressions  (3.11)  and  (3.15)  shows 
that  the  trapezoidal  rule  is  more  accurate  than  the  midpoint  rule  for  nonlinear 
problems  provided 

!ii  - 0 ( jy ' ) (3.16) 

The  two  one-step  formulas  have  been  used  to  integrate  the  model  nonlinear 
problem  (see  Park  [is] ) 

u + 100  tanh  (u)  = 0,  u(0)  = 0,  u(0)  = 25  (3.17 

The  results  show  that  the  midpoint  rule  gives  stable  solutions  for  all  stepsize 
ranges  tested  while  the  solutions  by  the  trapezoidal  rule  becomes  unstable 
for  stepsizes  h > 0.2.  The  interested  reader  is  referred  to  recent  work  by 
Danlquist  1 2 3 1 for  the  stability  analysis  of  multistep  one-leg  formulas. 

3 . 3 Step  Changing  Technique 

A proper  choice  of  implicit  integration  formulas  can  ensure  numerical  stability 
without  stepsize  limitations  (the  so-called  A-stability) . The  stepsize  is 
then  established  solely  from  accuracy  considerations.  Exploitation  of 

this  favorable  property  is  often  difficult,  however,  because  of  the  computa- 
tional overhead  associated  with  stepsize  changes.  These  difficulties  are 
examined  in  Section  3.3.1,  which  is  followed  by  the  presentation  of  an  improved 
matrix  scaling  technique  [24],  This  technique  accommodates  stepsize  changes 
without  requiring  a recalculation  of  the  matrix  Z in  (2.15)  . Applicable 
ranges  of  this  technique  and  the  importance  of  associated  predictors  are 
examined  in  Section  3.3.3.  Finally,  the  effect  of  stepsize  changes  on  the 
stability  of  the  integration  formula  is  discussed  in  Section  3.3.4. 


3-7 


3.3.1  Difficulties  in  Variable-Step  Implicit  Methods 


We  recall  Equations  (2.15)  and  (2.13)  for  the  reader's  convenience 


E u 
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The  stability  properties  of  the  integration  formula  (2.1')  can  be  maintained 
by  solving  (3.17a)  accurately  by  a convergent  iterative  method.  If  the  mass 
matrix  M is  diagonal,  the  following  predictor-corrector  iteration  schene 
appears  at  first  sight  attractive: 

u(m+1)  = M_1  [Ym)  - 6(D  + 5 K)  ,(m)l  (3.. 

— n *+  I -*n  ~ ~ — n J 

(This  results  from  splitting  JE*M  + (6  JO  + 5 2 K)  and  transferring  the  second 
term  to  the  right-hand  side.)  Unfortunately,  it  is  well  known  that  (3.18) 
converges  only  if 


cm  h < 1 

max 


(3.19) 


where  c is  a constant  of  order  unity.  The  stepsize  restriction  (3.19)  is 
of  the  same  order  as  that  imposed  by  explicit  methods.  Thus,  (3.18)  wipes 
out  the  basic  advantage  of  implicit  integration.  The  A-stability  of  implicit 
formulas  can  be  retained  by  using  the  Newton-like  iteration 


(ra)  , / (m+1) 

E (<5)  / u -• 

~ \ — n 


(3.20) 


Notice  that  for  linear  problems,  i.e.,  constant  E and  g (3.20)  gives  the 

^ (0) 

converged  solution  with  one  iteration  regardless  of  the  predictor  u 

n 

If  the  stepsize  is  varied,  the  matrix  E(5  ) has  to  be  refactored.  This 

~ new 

can  be  costly  in  practice  and  is  a major  barrier  to  the  effective  use  of 
implicit  formulas  for  variable  step  integration  of  the  equations  of  motion. 
A means  to  overcome  this  difficulty  is  described  below. 
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3.3.2  Matrix  Scaling  Technique 


Let  us  express  the  matrix  £(£)  in  the  font 


Ef<5)  = - E ( 5 ) - R 

— * s ~ o ~ 


(3.21) 


where  s is  a scaling  parameter  and  R is  a residual  matrix  to  be  determined. 
Substitution  of  (3.21)  into  (3.17a)  yields 


E (6  ) u - s(g  + Ru)  = C 

~ o — n -^n  ~ 


(3.22) 


where  the  residual  matrix  R is 


a ■ (j-  l)a  *(r  - 8)  s * (V- {!)  £ 


(3.23) 


An  iterative  method  for  (3.22)  takes  the  form 


(m+1) 
— n 


s E (5  )_1  (g  + R utm)  ) 
~ o ^n  — -n 


(3.24) 


Notice  that  we  have  freedom  for  improving  this  iterative  method  by  choosing 
>-he  scaling  parameter  s and  predictor  u/“‘J  appropriately.  The  convergence 
of  the  iterative  method  (3.24)  can  be  conveniently  examined  using 

the  scalar  case 


M = 1,  D = 2£w  and  K = w“ 


(3.25) 


The  spectral  radius  of  the  iteration  matrix  becomes 

| (1-s)  + 2(6  -s5)  + (62-s<52)  w‘ 


k [s  E_1(6  ) R]  = 


1 + 2 + u252 

o o 


(3.26) 


a)  Fixed-step  case  (s-1,  5^  - 6).  The  spectral  radius  is  zero  from  (3.26) 
and  accordingly  the  solution  converges  with  one  iteration  as  it  should. 


b)  Variable  Step  Case  (6^  f 6 ) . This  corresponds  to  integrating  with  a 

stepsize  other  than  the  factored  stepsize  5 . If  the  svstem  is  undamped 

o 

(5  = 0) , the  optimal  scaling  parameter  is 
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1 + 6 2 u:2 


1 t ? 2 
1 + 4-  u*- 


(3.27) 


F^r  the  multidimensional  case,  the  choice  (3.27)  with  u = u , a "scaling 

c 

frequency",  gives 


(H2  - 1)  (S3  - S32)  | 

C 

(1  + S32)  (1  + n 2V22 ) 


(3.28) 


in  which  n = 5 /5  , S3  = 6u. 

o 


The  case  S3  = 0 and  n 1 was  investigated  in  [25J  and  some  preliminary 
performance  data  was  reported  in  [5].  If  n < 1,  (3.23)  shows  that  ^ S2  ^ 
is  probably  adequate . 


Global  convergence  for  the  entire  frequency  range  is  obtained  with  the  following 
choice  of  the  scaling  parameter: 


(n2  - 1)  a2 

1 + rp  Q2 


(l  - n2) (fl2  - n2) 

max 

(1  + S32  ) (1  + n2S32) 

max 


s * 1 , n 51  1 


1 + 


1 + fl2 

max 


(3.29) 


, n < i 


As  can  be  seen  from  (3.29a),  the  spectral  radius  < for  n 1 rapidly  approaches 
unity  as  the  sampling  frequency  S3  increases.  This  restricts  the  choice  of 
predictors  to  extrapolation  from  displacements  only,  because  the  inclusion  of 
velocity  terms  in  the  predictor  can  generate  large  errors  in  the  predicted 
displacements  which  contain  high-frequency  components. 

On  the  other  hand,  when  n = 5/5  < 1,  the  spectral  radius  of  (3.29b)  approaches 

o 

unity  and  zero  for  low  and  high  frequency  components,  respectively.  A desir- 
able convergence  property  is  to  obtain  as  small  a spectral  radius  as  possible 
for  the  "accuracy"  range  0 < wh  < 2 and  to  bound  the  spectral  radius  below 
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unity  for  the  "noise”  range  wh  > 2 for  iteration  stability.  Note,  for 
example, that  the  globally  convergent  choices  (3.29)  yield 


(n2  - 1) 

1 + n2 


(l  - n2)  (-Q2  - l) 

max 


(i  + ) (i  + 

max 


r)  >>  1 


i 


(3.30) 


1 , n <<  1 ; SJ  ■+  1 


This  indicates  that  both  cases  give  a poor  convergence  rate  in  the  accuracy 

range  as  the  integrating  stepsize  <$  deviates  further  from  the  factored  stepsize 

5 . 

o 

3.3.3  Improvement  of  Low-Eh. . ~uency  Convergence  Rate 


In  order  to  improve  the  convergence  rate  for  the  accuracy  range  0 <_  oh  1 , 
we  introduce  an  acceleration  scheme 


(m+1) 


where 


« »<">  ♦(!-«)  u1"-11  - 8 r(n+l’ 


(m+1)  -1  , _ , . (m) 

L m 1 - s E ^ (2.+  R(s)  u ) 


(3.31) 

(3-32) 


Combining  (3.31)  with  (3.24)  one  can  derive  the  following  error  equation 


, (m+2)  (m+1)  , (m)  -1  (m+1) 

= a r + (1  - a r - 3 s E E Rr  (3.33) 

— — ~ ~o  ~ — 


The  stability  and  convergence  properties  of  the  acceleration  scheme  (3.31) 
can  be  evaluated  from  its  associated  characteristic  equation 

X2  I - (a  I - 6 s E E'1  R)  X+  (a  - 1)  I = 0 (3.34) 

It  can  be  shown  [24]  that  stability  requires  that 
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The  conditions  (3. 35) are  satisfied  by  the  choice 


U + U . 
max  min  „ 0 

a = t 6 , S = 


» mm  * max 


(3.36) 


where  y's  are  the  eigenvalues  of  the  system 


-1 


(s  E (6)  £(5  ) R - y I)  y * 0 


(3.37) 


The  convergence  rate  < of  the  acceleration  scheme  (3.31)  can  be  expressed  as 


k = 


1 - , 

Ju 

/u  . 

V max 

min 

1 + \fu  7u~ 

y max  mi 


(3.38) 


min 


Consequently,  the  narrower  the  spread  of  the  eigenvalue  spectrum,  the  faster 
the  convergence  rate. 


when  the  integration  stepsize  <$  is  larger  than  the  factored  stepsize  6 , i.e., 

o 

n = j /5  < 1,  one  obtains 
o 


(1  - n2)  a2 
It  a2 

(l  - n2)  p 

1 + en2 

0.17 (l-n2)/n2 


n ■+  o 


n -*•  ® 


Of  + 0.44/nz  ; n n = 1 
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(3.39) 


where  we  have  chosen  the  scaling  parameter  s and  the  constant  a as 


1 + n2a2 


1 + a 
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(3.40) 


(3.41) 

The  preceding  asymptotic  estimates  of  the  eigenvalues  of  (3.37)  indicate  that 

if  Q is  accurately  known  and  fi  6/6=1,  one  can  achieve  a rapidly  con-  ; 

max  max  o | 

vergent  scheme  through  the  acceleration  (3.31). 

A similarly  favorable  convergence  rate  for  the  case  n > 1 is  harder  to  achieve 

and  requires  alternating  iterations  with  two  different  factored  matrices  [24].  > 


= (1  + p n2)  n 


max 


3.34  Effect  of  Step  Changes  on  Stability 

There  are  two  techniques  for  incorporating  step  changes  in  the  multistep 


J 
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integration  formulas:  interpolation  of  past  solutions  to  represent  solutions 
at  equal  intervals  of  the  present  stepsize  and  use  of  variable  coefficients 
in  the  integration  formulas.  In  general,  interpolation  requires  additional 
solution  vectors  to  be  stored  for  preservation  of  the  fixed-step  accuracy  and 
stability.  The  introduction  of  variable  coefficients  alleviates  the  addi- 
tional storage  requirement,  but  is  known  to  cause  instability  if  the  stepsize 
is  varied  too  often  [26] . Figure  6 shows  the  stability  region  of  a three-step 
formula  [l8]  as  adopted  to  variable  coefficient  technique  in  the  form 


where 


a u . = a u + a_  u . + a,  u _ + h u 

o n+1  In  2 n-1  3 n-2  n+1 
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1 + n. 


2 (1  + n,  + n2) 
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2n, 


(nx  + n2) 
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(1  + 3^) 
(l  + n1) 
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(1  + T^J 


(l  + n + n2)  (nL  + n2)  n2 


n.  = h /h  , = h /h  . . 

1 n n+1  2 n-1  n+1 


(3.42) 


Note  from  Figure  6 that  a decrease  in  the  stepsize  shrinks  the  unstable  zone 
and  results  in  an  increased  numerical  damping  as  the  stability  boundary  departs 
further  from  the  imaginary  axis.  An  increase  in  the  stepsize,  on  the  other 
hand,  expands  the  unstable  zone  which  includes  a segment  of  the  imaginary  axis. 
This  segment  will  cause  local  instability  if  the  combination  of  frequency 
contents  and  the  stepsize, 
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Fig.  6 Effect  of  Uneven  Step  Interval  on  Local  Stability 


3-15 


oh,  falls  in  that  zone.  It  is  also  noted  that  the  fixed,  k-step  formula  is 
-ecovered  at  k-th  step  after  the  stepsize  is  changed. 

How  seriously  stepsize  changes  affect  the  global  stability  and  accuracy  has 
not  been  systematically  studied  yet.  Experience  indicates  that,  for  k-step 

formulas,  the  effect  is  insignificant  if  the  stepsize  is  varied  at  most 
every  2 k-th  step. 


3.3.5  Step  Changing  Techniques  for  the  Central  Difference  Formula 

The  effect  of  step  changing  techniques  on  the  stability  of  explicit  formulas 
is  more  pronounced  than  on  implicit  formulas  because  the  stability  margin  of 
explicit  formulas  is  quite  narrow,  mostly  wh  <_  2.0.  We  present  a relevant 
portion  of  the  results  from  [27]  for  the  central  difference  formula.  Of 
the  four  techniques  considered  (see  Table  5) , the  velocity  averaging  technique 
(FAVE)  is  most  desirable  for  undamped  problems,  whereas  the  variable  coeffi- 
cient technique  (FVCO)  outperforms  the  rest  for  heavily  damped  problems.  As 
in  the  case  of  implicit  formulas,  the  question  of  how  often  the  stepsize  can 
be  changed  without  a noticeable  affect  on  the  global  solution  accuracy  has 
not  been  satisfactorily  resolved  yet. 

3 . 4 Stepsize  Selection  Strategies 

The  earliest  known  numerical  integration  method  is  due  to  Euler  [23,29]. 

By  1973,  the  number  of  publications  had  grown  to  over  600  [30]  and  over  50 
articles  are  published  every  year.  Of  these,  at  least  a dozen  are  directly 
concerned  with  the  class  of  stiff-oscillatory  systems  to  which  structural 
dynamic  problems  belong.  Adequate  stepsize  selection  strategies  for  non- 
stiff parabolic  problems  exist  [3l]  and  some  promising  strategies  for  stiff 
parabolic  problems  have  been  published  [32].  The  authors  are  not  aware, 
however,  of  many  publications  that  deal  with  stepsize  selection  for  stiff 
oscillatory  problems. 

The  paucity  of  publications  is  apparently  due  to  the  fact  that  low-accuracy 
solutions  (e.g.,  two  correct  digits)  are  acceptable  for  many  structural 
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Table  5 


FORMULAS  FOR  STEP  CHANGES 


3-17 


’ '■'$ 

'it 

j.% 


dynamic  problems,  while  most  of  the  available  stepsize  selection  strategies 
are  based  on  the  truncation  error  concept,  which  is  effective  chiefly  for 
high-accuracy  requirements. 

For  variable-step  explicit  integration,  the  ability  of  detecting  local 
(incipient)  instability  is  mandatory.  For  implicit  integration,  the  ability 
of  maintaining  a given  accuracy  level  is  essential.  A general-purpose  inte- 
gration package  should  therefore  have  the  ability  to  monitor  both  cases 
depending  on  the  integration  mode,  and  transition  from  one  mode  to  the  other 
should  be  governed  by  anticipated  cost  minimization.  We  now  examine  these 
goals  in  more  detail . 

3.4.1  Explicit  Variable-Step  Strategies 

In  low— accuracy  explicit  integration,  the  detection  of  instability  becomes 

a major  concern.  This  requires  an  accurate  assessment  of  the  highest  fre- 
quency of  the  system.  Experience  shows  that  this  is  not  a simple  task  in 
nonlinear  problems,  and  repeated  calculations  to  that  effect  can  be  costly. 
Numerical  experiments  have  shown  that  monitoring  the  truncation  error  as 
stability  guard  is  generally  useless.  Of  several  empirical  techniques 
proposed  and  tested  so  far,  the  following  "perturbed  apparent  local  frequency" 
concept  has  been  found  to  be  reliable  and  efficient.  The  perturbed  apparent 
frequency  is  defined  as 


Au . Au . 
i 1 


(3.43) 


which  can  be  expressed  for  linear  homogeneous  equations  as  [27] 
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where  t_.  is  the  i-th  column  eigenvector  and  £ -s  the  uncoupled  modal  displace- 
ment vector  from  the  relationship 


u 


(3.45) 
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Remark  1 . The  qualifier  "apparent"  is  used  here  to  emphasize  only  those 
frequencies  that  are  participating  in  the  dynamic  response.  For  example, 
a structure  under  free-fall  gravity  motion  exhibits  only  the  rigid-body 
motion  (zero  frequency).  In  this  case,  is  zero  and  any  arbitrarily 

large  stepsize  shculd  yield  an  exact  solution,  which  is  confirmed  by  numeri- 
cal experiments  [33]. 

Remark  2 . The  highest  frequencies  participate  in  the  dynamic  response  at 
early  time  and  may  hardly  participate  at  a later  time.  The  apparent  fre- 
quency concept  should  in  this  case  result  in  an  increased  stepsize  at  the 
late  time  period,  which  appears  also  to  be  the  case. 


Numerical  experiments  indicate  that  the  step  selection  strategy  based  on  the 
apparent-frequency  sampling  theory  outperforms  all  other  strategies 

tested.  Further  refinements  cf  this  strategy  are  expected  by  experimenting 
on  large-scale  problems. 

3.4.1  Implicit  Variable-Step  Strategies 

Four  strategies  for  stepsize  selection  in  implicit  integration  will  be 
mentioned.  They  are  based  on  the  local  truncation  error,  residual  term, 
Richardson's  extrapolation  and  deviation  from  linearity,  respecu  vely. 


The  use  of  the  truncation  error  is  well  known  and  need  not  be  further 
explained,  see  e.g.  [7].  The  residual  term  d is  the  difference  between  the 
acceleration  computed  from  the  differential  equations  of  motion  (1.1)  and 
that  computed  from  the  integration  formula  (2.7b): 


d = u(u  ,u  ) - (u  - h )/h. 
— n n — n ~n  ~n  5 


(3.46) 


An  approximate  local  error  e can  be  expressed  [34]  as 
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where  p is  the  order  of  the  integration  formula. 


The  strategy  based  on  Richardson's  extrapolation  uses  the  difference 
between  two  solutions  obtained  with  stepsizes  h and  y h,  respectively,  to 
express  the  local  error  as 

e ^ u (r  h)  * a (h)  (3.48) 

-n  — n z -n 

or  a similar  estimate  of  when  smoothing  is  performed  [33]  . Limited 
experience  suggests  that  the  residual  term  estimate  (3.47)  is  far  more  reliable 
than  the  truncation  error.  The  estimate  (3.48)  is  the  most  reliable  of  the 
three,  but  requires  at  least  three  solutions  per  step.  In  future  numerical 
experiments,  error  estimation  by  extrapolating  at  every  tenth  step  is  con- 
templated to  reduce  the  computational  overhead . 

Finally,  the  deviation  from  linearity  can  be  used  as  stepsize  selection 
strategy.  This  works  effectively  if  the  major  source  of  error  is  due  to 
system  nonlinearities  and  one  has  a priori  knowledge  of  the  linear  behavior 
of  the  system.  This  option  is  attractive  in  those  cases  where  the  analyst 
performs  a linear  dynamic  analysis  first  and  then  decides  to  proceed  to  a 
nonlinear  analysis  on  account  of  critical  design  considerations. 

No  consensus  as  to  which  strategy  is  most  suitable  for  structural  dynamic 
analysis  has  been  reached  as  yet.  In  fact,  this  constitutes  an  important 
area  of  near-future  research  (cf.  Section  3.6). 

3 . 5 New  Approaches 

It  should  be  clear  from  previous  sections  that  the  relative  performance  of 
explicit  and  implicit  methods  depends  on  problem  characteristics,  user 
requirements  (especially  accuracy  level) , and  the  computer  cost  function. 

The  two  methods  can  be  applied  concurrently  over  different  portions  of  the 
structure  [36,37]  or  in  an  interleaving  manner  so  that  the  integration  pro- 
cess alternates  between  the  two  modes  according  to  a cost  minimization 
criterion  [38], 
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A combination  of  explicit  integration  and  mode  superposition  methods  can 
be  used  to  increase  the  effective  integration  stepsize.  Knowledge  of  a set 
of  high  frequency  dominant  eigenvectors  allow  the  corresponding  subspace  to 
be  annihilated  from  the  equations  of  motion  [39].  The  net  result  of  this 
technique  is  an  extension  of  the  maximum  stable  stepsize  for  explicit  methods; 
this  gain  may  be  significant  for  certain  classes  of  problems. 

It  may  also  be  possible  to  combine  the  traditional  Rayleiah-Ritz  method 
with  implicit  methods.  Direct  time  integration  would  be  used  only  occasionally 
to  evaluate  errors  in  the  Rayleigh-Ritz  approximations  and,  if  necessary, 
to  trigger  the  injection  of  additional  sets  of  trial  mode  shapes.  Such  an 
effort  would  probably  shed  some  light  on  the  usefulness  of  the  so-called 
exponentially  fitted  formulas  from  an  engineering  application  viewpoint. 

The  classes  of  problems  for  which  each  of  the  aforementioned  approaches 
become  advantageous  has  not  been  fully  explored,  and  more  extensive  testing 
will  be  required  before  definite  conclusions  to  this  effect  are  reached. 

3 . 6 Areas  for  Future  Research 


The  increasing  importance  of  dynamic  analysis  in  design  and  verification  of 
complex  structures  has  given  renewed  impetus  to  research  activities  in 
direct  time  integration  methods.  During  the  past  decade,  significant  progress 
has  been  made  in  the  following  areas:  techniques  for  stability/accuracy 
assessment  in  linear  problems,  qualitative  assessment  of  effect  of  nonlineari- 
ties on  stability,  comparative  performance  rating  of  multistep  formulas  on 
various  nonlinear  model  problems,  and  spatial  hybridization  of  explicit  and 
implicit  formulas. 

Three  important  near-term  research  areas  will  be  mentioned.  The  authors 
believe  that  the  most  important  one  is  the  further  development  of  stepsize 
control  strategies  for  implicit  integration  methods.  Progress  in  this  area 
would  facilitate  the  production  of  cost  effective  dynamic  analyzers  that  can 
operate  with  the  same  degree  of  reliability  as  that  of  presently  available 
static  analyzers,  and  hence  attract  interest  from  practicing  engineers. 
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A second  area  involves  the  systematic  study  of  performance  crossover  points 
between  explicit  and  implicit  methods  for  a given  computer  environment  and 
a set  of  representative  problems.  The  crossover  pcint  can  be  expected  to 
depend  not  only  upon  the  required  accuracy  level  and  the  stepsize  selection 
technique,  but  also  upon  the  computer  characteristics  and  the  quality  of  the 
supporting  software.  Resolution  of  this  topic  would  simplify  the  implemen- 
tation of  problem- adaptive  temporal  switches  between  both  integration  modes. 

A third  area  pertains  to  the  study  of  the  influence  of  formulations  and 
computational  procedures  for  nonlinear  problems  on  the  stability  of  the 
integration  formulas.  (In  other  words,  investigation  of  the  interaction 
between  Sections  2 and  3.2.) 

Finally,  as  a long-term  research  topic,  we  should  mention  the  potential 
impact  of  new  computing  equipment  on  the  performance  of  direct  time  inte- 
gration methods . Some  of  the  techniques  which  are  currently  deemed  obsolete 
may  turn  out  to  be  most  suitable,  or  new  formulas  and  procedures  may  have  yet 
to  be  developed  to  fit  new  computing  environments.  At  the  present,  however, 
the  development  of  new  integration  formulas  is  not  likely  to  advance  the 
state  of  the  art  as  significantly  as  that  of  the  three  near-term  areas. 


4.  CONCLUDING  REMARKS 


Substantial  computational  cost  reductions,  of  the  order  of  2 to  50, have  been 
observed  as  a result  of  the  use  of  the  generalized  pseudo-force  approach  over 
the  conventional  tangent  stiffness  approach.  Another  cost  reduction  of 
roughly  2 to  5 can  be  achieved  by  paying  careful  attention  to  implementation 
details.  The  total  gain  can  thus  span  one  or  two  orders  of  magnitude. 

An  important  consequence  is  that  nonlinear  dynamic  analysis  can  no  longer  be 
regarded  as  being  orders  of  magnitude  more  expensive  than  the  corresponding 
linear  analysis.  For  sufficiently  large  problems  the  array-processing  effort 
involved  in  the  solution-advancing  cycle  dominates  over  that  of  evaluating 
nonlinear  terms.  Under  such  circumstances,  and  assuming  identical  inte- 
gration-driving software,  the  (nonlinear/linear)  cost  ratio  per  time  step 
is  close  to  unity  in  explicit  integration.  For  implicit  integration,  the 
ratio  depends  primarily  on  the  number  of  iteration  cycles  per  step  if  refac- 
torizations are  avoided,  and  typically  varies  from  2 to  5.  These  significant 
cost  reductions  have  been  translated  into  Figure  7;  this  is  an  update  of  a 
chart  prepared  in  1976  [3] , which  reflected  the  then  prevalent  use  of  the 
tangent  stiffness  approach. 

It  is  hoped  that  the  improved  computational  efficiency,  when  coupled  to  the 
increased  reliability  brought  about  by  the  problem-adaptive  techniques 
discussed  in  Section  3,  will  eventually  result  in  breeder  usage  of 
transient  nonlinear  analysis  by  practicing  engineers 
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Appendix  A 

COMPUTER  IMPLEMENTATION  OF  DYNAMIC  ANALYSIS 


A. I Organizational  Approaches 

Two  philosophies  may  be  followed  in  the  implementation  of  transient  response 
analysis  capabilities. 

1.  Tightly  Coupled  Processors.  The  analysis  control  is  embedded 
within  the  framework  of  a specific  structural  anal/zer. 

2.  Loosely  Coupled  Processors.  The  analysis  control  is  segregated 
as  a modular  "driver"  processor,  which  is  interfaced  with  existing 
analyzers  through  a data  base  management  system. 

(A  similar  choice  appears  in  other  application  programming  areas  such  as 
optimization  and  synthesis,  in  which  a general  purpose  driver  processor 
controls  analysis  processors.) 

The  first  approach  offers  the  benefits  of  expedience  and  potentially  high 
computational  efficiency,  inasmuch  as  special  analyzer  features  can  be 
exploited.  It  is  indicated  if  the  development  is  carried  out  by  the 
designer  of  the  analyzer,  or  a tightly  knit  team  that  includes  the 
designer.  The  main  disadvantage  of  this  approach  is  the  sacrifice  in 
flexibility  of  application,  main '.ainability  (the  ability  to  repair  errors)  , 
and  modifiability  (the  ability  to  make  functional  changes) . Furthermore, 
if  the  structural  analyzer  is  * ery  complex,  or  is  adjudged  proprietary  by 
its  source,  this  approach  may  be  unfeasible. 

The  second  approach  is  consistent  with  present  trends  as  regards 
construction  of  data  base  linked  program  networks  [40] . Consequently,  it 
is  believed  to  offer  the  greatest  potential  on  a long  term  basis.  It  is 
particularly  appealing  when  several  existing  analyzers,  some  of  which  may 
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be  unmodif iablo , are  to  be  interfaced  with  the  same  integration  driver. 

One  drawback  of  this  approach  is  the  need  for  considerable  attention  on 
the  part  of  the  developer  to  the  overall  design  of  the  program  to  avoid 
"downstream"  surprises.  Separable  driver  operation  also  presupposes  the  ready 
availability  of  a library  of  supporting  matrix  processors  as  well  as  that  of 
a scientific  data  management  system.  Because  of  the  long-term  nature  of  the 
development  effort,  this  approach,  which  i9  the  one  discussed  in  the  following 
sections,  is  not  necessarily  suitable  for  every  programming  group. 

A. 2 The  Plavers 


Many  industrial  and  engineering  organizations  maintain  a library  of  structural 
analyzers.  These  are  either  acquired  from  external  sources  or,  if  the 
organization  is  sufficiently  large,  developed  in  house.  Most  of  the  analyzers 
will  be  linear,  some  nonlinear;  in  the  latter  case,  a restricted  problem 
application  range  is  the  norm  rather  than  the  exception. 

Another  intervening  software  element  is  the  time  integration  control  module, 
herein  called  the  integrator  (I)  for  brevity. 

The  function  of  the  analyzer  (A)  is  to  generate  the  discrete  dynamical 
equations,  i.e.,  the  components  of  (1.1).  The  integrator  advances  the 
state  solution  by  following  one  of  the  computational  sequences  discussed 
in  Section  2.  Computed  responses  may  be  fed  back  to  the  analyzer  to  regene- 
rate nonlinear  terms  or  to  compute  derived  quantities  such  as  stresses. 

Upon  analysis  completion,  the  computed  responses  are  sent  to  a display  post- 
processor (D) . This  is  usually  a separate  module,  although  portions  of  D 
may  be  actually  duplicated  in  the  integrator  to  produce  "display  snapshots" 
during  the  analysis  process. 

The  process  just  described  is  schematized  in  Figure  A-l.  This  software 
organization  may  be  regarded  as  an  instance  of  modularization  by  program 
function  [ 4l] - 

A. 3 Transient  Linear  Analysis 


A large  percentage  of  dynamic  analysis  concern  linear  problems.  Moreover, 


Fig.  A.rl  Basic  Software  Elements  in  Structural 
Dynamic  Analysis 


an  important  class  of  nonlinear  problems  can  be  processed  in  "piecewise 
linear"  mode,  in  which  the  response  is  linear  except  at  isolated  "busy" 
periods.  Consequently,  the  development  of  a modular,  general-purpose 
integrator  would  be  perhaps  economically  unjustifiable  if  a capability  for 
efficiently  functioning  in  linear  or  piecewise-linear  mode  is  not  provided. 
This  capability  is  activated  by  suppressing  nonlinear  feedback. 


Figure  A-2  schematizes  the  operational  flow  for  transient  linear  analysis. 

In  this  case,  a high  degree  of  modularity  is  passible.  The  structural 
analyzer  functions  primarily  as  an  array  generator.  These  arrays  typically 
include  the  matrices  K,  D,  and  M and  the  "base  load  vectors"  L.  discussed 

. ~ ~ ■—  -i 

below.  These  arrays  are  placed  on  a global  data  base  (GDB) , which  resides 
on  permanent  storage  devices.  If  the  formatting  requirements  imposed  by  the 
support  processing  utilities  (F)  conflict  with  those  of  the  analyzer,  the 
arrays  must  be  passed  through  a data  format  converter  (C^)  before  they  are 
ready  to  be  submitted  to  the  integrator. 

The  temporal  variation  of  the  applied  force  vector  can  be  often  expressed 
in  the  separable  form 

f(t)  = ^ li(t)  = L \ (A.l) 

i=l 

where  are  user  supplied  load  factors  and  L is  the  base  load  matrix  whose 
columns  are  the  base  load  vectors  L. . Frequently,  n.  <<  n , the  number  of 
freedoms;  if  = 1,  the  loading  is  called  proportional.  Usually,  the  base 
load  vectors  result  from  discretization  of  forcing  fields  such  as  surface 
tractions,  temperature  variations  or  body  forces,  and  must  be  therefore 
generated  by  the  structural  analyzer  (this  is  the  cause  of  the  path  accessing 
the  applied  force  segment  in  Figure  A-2).  The  decomposition  (A.l)  entails 
considerable  computational  advantages  if  L is  constant  over  a set  of  problems 
and  che  user  wishes  to  investigate  the  response  to  a ensemble  of  load 
histories . 


computed  responses  are  placed  on  the  global  data  base  and  can  be  later 
accessed  by  the  display  postprocessor  module  (D) . If  the  history  of  derived 
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quantities  such  as  strains  and  stresses  is  required,  the  simplest  procedure 
is  to  send  the  computed  responses  back  to  the  analyzer;  derived  responses 
can  be  then  accessed  by  D.  Two  format  converters,  identified  as  and 
in  Figure  A-2,  are  generally  required  for  this  operation.  Displayed  results 
may  take  the  form  of  time  histories,  envelope  (peak  response)  plots,  and 
response  spectra. 

An  implementation  of  the  operational  flowchart  in  the  form  of  data  base-linked 
program  network  [4o]  is  illustrated  in  Figure  A-3.  A key  ingredient  is  a 
common  data  base  manager  (M)  through  which  the  global  data  base  is  accessed. 

In  addition,  computational  modules  such  as  the  analyzer  and  the  supporting 
matrix  processors  normally  have  their  own  local  data  bases  (LDB) , which  con- 
tain intermediate  result  and  auxiliary  data  structures.  These  local  data 
bases  are  efficiently  maintained  in  a "virtual  memory"  storage  pool  adminis- 
tered through  local  data  base  managers.  These  may  be  developer-supplied  or 
part  of  the  operating  system  library. 

It  is  important  to  realize  that  the  degree  of  modularity  of  this  implementa- 
tion is  very  high.  For  example,  replacing  an  analyzer  by  another  entails 
changing  the  converters  (C.  , C,,  C ),  with  no  side  effects  propagating  to 
the  integrator.  Conversely,  incorporation  of  technical  improvements  in 
the  integrator  does  not  imply  that  the  analyzers  must  be  modified. 

Th  ability  to  "plug  in"  different  linear  analyzers  is  particularly  advanta- 
geous if  segments  of  the  structure  are  modeled  by  different  programs.  This 
is  a common  occurrence  in  large-scale  engineering  projects  involving  many 
subcontractors.  The  converter  is  then  complemented  by  a "substructure- 
linking"  module  that  ties  the  various  models  together. 

Another  advantage  of  the  data  base-linked  organization  is  potentially  smooth 
adaptability  to  distributed  processing  (segments  of  the  problem  may  be 
assigned  to  different  computers) . For  example,  the  display  postprocessing 
stage  can  be  most  economically  carried  out  on  a minicomputer. 

A. 4 Transient  Nonlinear  Analysis 

Figure  A-5  depicts  the  operational  flow  for  transient  nonlinear  analysis 
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carried  out  through  a pseudo-force  based  integrator.  The  key  difference 
froin  the  linear  case  is  the  occurrence  of  response  feedback  at  each  time 
step.  This  has  a detrimental  effect  on  the  degree  of  modularity  attainable 
as  regards  integrator-analyzer  interaction. 

Computed  response  information  is  used  by  the  analyzer  for  evaluating  stresses, 
nonlinear  force  terms  and  applied  forces  (if  the  latter  are  configuration- 
dependent)  . In  problems  involving  material  nonlinearities,  the  stress  cal- 
culations also  impacts  the  generation  of  stiffness  matrices  and  of  certain 
nonlinear  force  terms . 

One  of  the  constraints  imposed  by  the  lessened  functional  modularity  is 
the  elimination  of  the  array  format  converter  C . The  array  conversion 
would  have  to  be  exercised  many  times  during  the  run,  which  can  be  expensive. 
The  elimination  of  can  be  effected  in  several  ways: 

1.  Embed  in  A so  that  arrays  are  directly  produced  in  P-acceptable 
format. 

2.  Configure  P so  that  various  self-labeled  array  structures  are 
accepted  [42], 

3.  Do  matrix  processing  in  A so  that  most  of  the  required  array  trans- 
mission is  in  the  form  of  one-dimensional  arrays. 

In  connection  with  the  third  approach,  it  is  noteworthy  to  mention  that 
explicit  nonlinear  integrators  are  inherently  more  modular  than  implicit 
integrators,  inasmuch  as  all  analyzer-to-integrator  data  can  be  transported 
in  vector  form,  including  the  diagonal  mass  matrix  M (cf.  Section  2.7). 

How  can  the  functional  flowchart  of  Figure  A-4  be  implemented  to  maximize 
modularity?  One  possible  solution  is  shown  in  Figure  A-5.  Here  the  analyzer 
and  integrator  operate  as  "coroutines"  linked  by  an  executive  control 
component  E. 

Some  operating  systems  require  that  E be  a developer-supplied  program,  a 
restriction  that  forces  A and  I-P  to  be  organized  as  parallel  overlays. 
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complete  reconstruction  of  the  overlev  image.  V.cre  auvar.ced  operating 
systems  allow  Z to  he  implemented  as  a catalogued  control  procedure  that 
simply  executes  A cr  1 according  to  run  condition  state  variables.  In 
this  case,  the  two  sides  of  Figure  A-5  can  he  kept  as  physically  independent 
programs  served  by  a common  version  of  the  data  base  manager 

A. 5 The  STINT  Code 


A stand-alone  integrator  called  STINT  has  been  developed  following  the  loosely- 
coupled  processor  approach.  Initial  versions  were  primarily  used  as  testbed 
of  new  techniques  and  reference  source  for  subsequent  development  of  produc- 
tion- level  versions.  STINT  car.  operate  as  an  implicit  or  explicit  integrator. 

The  implicit  integration  branch  uses  the  trapezoidal  rule  and  the  Park 
3-step  method  [l3]  as  basic  formulas  for  treating  linear  and  nonlinear 
problems,  respectively.  Processing  of  the  latter  relies  upon  a general 
pseudo-force  approach  used  in  conjunction  with  the  matrix  scaling  technique 
described  in  Section  3.  The  (Cl)  computational  sequence  is  presently  used, 
although  implementation  of  (JO)  is  anticipated.  Further  implementation  details 
are  given  in  Underwood  and  Park  [5]. 

The  explicit  integration  branch  is  based  upon  the  variable-step  central  dif- 
ference scheme  [27] . Implementation  details  are  given  in  Underwood  and  Park [3 3 ] 

A long-term  effort  in  "tuning  up"  problem-adaptive  features  of  STINT, 
notably  the  automatic  stepsice  selection,  has  resulted  in  impressive  but 
"staggered"  performance  gains  in  either  the  implicit  or  explicit  branch.  As 
a consequence,  efficiency  crossover  points  on  the  same  test  problem  set  have 
varied  as  new  verions  are  prepared.  Presently,  the  explicit  integration  branch 
is  considered  to  be  in  superior  shape. 

Although  STINT  can  switch  from  explicit  to  implicit  mode,  or  vice-versa, 
during  an  analysis  run,  the  switch  points  must  be  presently  specified  by 
the  user.  Cur  ultimate  goal  is  to  have  STINT  select  the  integration  mode 
automatically  according  to  a cost  (merit!  function  minimization  criterion,- 
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this  feature  would  greatly  enhance  the  problem  adaptivity  of  the  integrator. 
The  realization  of  this  goal  will  have  to  wait,  however,  until  stepsize 
control  techniques  in  the  implicit  branch  demonstrate  the  same  degree  of 
reliability  currently  achieved  in  the  explicit  branch. 
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